A cluster algorithm for resistively shunted Josephson junctions 
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We present a cluster algorithm for resistively shunted Josephson junctions or similar physical 
systems, which dramatically improves sampling efficiency. The algorithm combines local updates 
in Fourier space with rejection-free cluster updates which exploit the symmetries of the Josephson 
coupling energy. As an application, we consider the superconductor-to-insulator transition in a 
single junction at intermediate Josephson coupling and determine the temperature dependence of 
the zero bias resistance as a function of dissipation strength. 



Resistively shunted Josephson junctions (JJ) at zero 
temperature undergo a superconductor-to-insulator tran- 
sition if the shunt resistance is increased beyond a 
critical value, which equals the quantum of resistance 
Rq — h/Ae^. This dissipative phase transition was first 
predicted by Chakravarty 1] and Schmid jSj and sub- 
sequently studied by several authors |^ |M IS - Simi- 
lar transitions occur in JJ arrays and rich phase dia- 
grams have been worked out using variational approx- 
imations [3, mean field theory and semi-classical cal- 
culations |9l| . Recently, experimental realizations of re- 
sistively shunted single junctions 10] and JJ-chains 
have allowed to test some of the theoretical predictions. 
These experimental advances and the rising interest in 
questions related to dissipation in macroscopic quantum 
devices have led to new theoretical investigations. Float- 
ing phases in linear arrays [l^ and complicated phase di- 
agrams for a two-junction system with charge relaxation 
|13| have been pred icted. An extension of the latter work 
to chains of JJ |l4| attempts to explain the intriguing ex- 
perimental results on the superconductor-to-normal tran- 
sition in nanowires . 

It would be desirable to test these sometimes contro- 
versial theoretical predictions numerically. Path-integral 
Monte Carlo studies of a single, resistively shunted JJ 
have recently been carried out jl^ fl7l |. However, the 
lack of efficient algorithms has so far proven to be an 
obstacle and the results in Ref. even disagree with 
analytical predictions. In this Letter, we present a new 
cluster algorithm which dramatically improves sampling 
efficiency compared to local update schemes and allows 
to simulate the junction at more than 100 times lower 
temperature. Such efficient updates for single junctions 
will be important for the simulation of JJ chains and 
arrays. 

Dissipation in the shunt resistor is introduced by cou- 
pling the phase difference (p across the junction linearly to 
an infinite set of harmonic oscillators [la . For an Ohmic 
heat bath one finds the effective action [6| (setting h—1) 
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The first term describes the capacitive coupling be- 
tween the two superconductors (the effective charging 
energy of the junction, Ec = sets the energy scale) 
and the Josephson energy for a coupfing strength Ej. 
The last term results from integrating out the heat bath 
degrees of freedom and introduces dissipation. Because 
of the Josephson energy and the non-compact nature of 
the phase variable one cannot directly employ the cluster 
algorithms available for spin systems ,2Qe , as was 
done in Ref. jl^ • We therefore propose a new algorithm 
consisting of two kinds of updates: (i) local updates in 
Fourier space compatible with the Gaussian terms in 
and (ii) rejection-free cluster updates. The fist type of 
moves assures ergodicity of the algorithm and the second 
type produces global cluster updates compatible with the 
energetic constraints from the Josephson term. 

We discretize imaginary time into N (assumed odd) 
time steps At — (3/N and introduce the dimensionlcss 
coefficient a — Rq/Rs- The action |^ can then be ex- 
pressed in the simple form 
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S[(j>] = ^ ak\(f)kf 
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EjAt J2 cos(0,), (2) 
j=o 

where (pk = X^jlJo^ e^^^^(j)j denotes the Fourier trans- 
form of (j). The positive coefficients ak are defined as 
o-k — jjign " 9k), with gk the Fourier transform of the 
kernel (j 7^ 0) 
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Since (j)k — (t>N-k, only {(f>k\k — 0, . . . {N ~ l)/2} need to 
be considered. In a local update of the frequency com- 
ponent (pk , we choose a new value according to the prob- 
ability distribution of the Gaussian term in Eq. ij^J, 



(4) 



and accept it with probability 

p(</'oid ^ 0„ew) = min (1, e-{S'^[*— . (5) 
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Here Sj[(t)] = -EjAtY^'^'^^^ cos{(j)j) and 4>oid, (t'new de- 
note the backward Fourier transform of the old and new 
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fc-space configuration, respectively. Such local updates 
can be performed in a time 0{N) and have recently been 
used in the simulation of 2D J J arrays [2]| . 

For reasonably large values of Ej, local fc-space up- 
dates which introduce phase changes on the order of 
2tt will be strongly suppressed, because their sinusoidal 
shape does not resemble an optimal phase slip path. Al- 
gorithms based on local updates alone will therefore be 
ineffective near the phase transition, where phase slips 
start to proliferate. A typical path will stay most of the 
time near one or the other of the minima in the cosine- 
potential, as shown in Fig. 1. Global moves compatible 
with this step- like structure would be updates which shift 
the phases 4>j by ±2tt in some random interval or reflec- 
tions about (f> — riTT, n G Z, that is the positions of the 
maxima and minima of the cosine potential. 

We exploit the latter symmetry to design a rejection- 
free cluster update consisting of the following four steps, 
which are illustrated in Fig. (i) An axis (f> = n^^^^ir 
with integral n^^^^ in the range [— rimax, J^max] is randomly 
chosen (the significance of nmax is discussed below) and a 
random site jroot is picked as the root site of the cluster. 
We introduce relative coordinates 

^f' ^ ^1 - (6) 

which are updated in a cluster move as (j)'^^^^ — > — 
Such updates do not change the value of Sj. (ii) The 
cluster of sites connected to jroot is constructed in a way 
analogous to other cluster algorithms 0, E^. Esl |26| . 
although the task is complicated by the presence of long 
ranged interactions. Two sites at positions i and j are 
connected with probability 

Pit, j) = max (0, 1 - e-{^[^r^-0r=l-S[*r".0r=l}) , (7) 
where 

(iii) The phases of the sites j belonging to the cluster 
are updated according to 0j 2n^^^^TT — (pj (in relative 
coordinates (p^^^^ —cf)^^^^). (iv) If necessary, the whole 
configuration is shifted by multiples of ±27r in such a way 
that the mean value (j> — 'Yl!j=o 't'j closest to the 
potential minimum corresponding to n'^^^^ = 0. The last 
step is important to assure detailed balance with a finite 
Hmax- The parameter rimax must be large enough that 
the re-centered configuration is contained in the interval 
[— 'T-maxTT, nmaxTr]. This valuc Can be determined during 
the thermalization of the system. An alternative to fixing 
the interval of symmetry points and re-centering the new 
configurations would be to choose the axis among the 
n-max symmetry points above or below ■ 

Using ideas of Luijten and Blote [l^, the cluster move 
can be performed in a time 0{N \ogN) despite long 
ranged interactions, which allows us to compute precise 
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FIG. 1: Illustration of the cluster algorithm for non-compact 
phase variables. The figures show from top to bottom: (i) 
Original phase configuration. Possible axis positions are indi- 
cated with dashed lines, located at the maxima and minima 
of the cosine potential. The randomly chosen axis and root 
site of the cluster are marked with the black solid line and 
black dot respectively, (ii) Cluster of connected sites. The 
sites which could potentially connect to the root site are indi- 
cated with tic marks, (iii) New phase configuration obtained 
by flipping the cluster around the axis. 

data for large systems (up to TV « 10'* if AtEc = 0.25). 
On small lattices the data produced by local update 
schemes are consistent with those obtained using the clus- 
ter algorithm, but we find that the results from local 
update simulations become unreliable for systems larger 
than N w 100 because of diverging autocorrelation times. 
Fig. 121 plots the integrated autocorrelation time r for 
{{4> — (f>)^) as a function of system size N. Even though 
the CPU time for a cluster update is 0{N log N) as com- 
pared to 0{N) for a local update, the gain in sampling 
efficiency is considerable. 

To demonstrate the performance of the new algorithm 
we study the localization transition at intermediate val- 
ues of the Josephson coupling energy. Theoretically, this 
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FIG. 2: Integrated autocorrelation times r for {{(f> — (j))^) ver- 
sus system size A^. The data were obtained at the critical 
point a — 1 using Ej/Ec = 1. 



transition is predicted to occur at a = 1 in the limit of 
large or small Ej/Ec and it was conjectured that its 
position is in fact independent of Ej/Ec- Local up- 
date Monte Carlo simulations have recently been used 
in Ref. jl6l | to investigate the transition at intermediate 
coupHng energy. These authors plotted ((0 — 0)^) as a 
function of dissipation strength a and claimed that the 
slope of such curves changes abruptly at the critical value 
Qfc = 1. In the inset of Fig. |3|we show our high accuracy 
data for the same parameter values as in Ref. 0| . Even 
at considerably lower temperatures, {{(j) — ^)^) changes 
smoothly around a ~ 1, and the method of Ref. [ig 
cannot be used to determine the critical point. 

Instead, we plot {{(j) — 4>)'^) as a function of the inverse 
temperature f3 (or system size N). For a — 1 one ob- 
serves an increase proportional to log A''. Below a = 1, 
the fluctuations grow faster than logarithmically and 
hence diverge (insulating phase), whereas above a — 1 
the phase fluctuations presumably saturate at some finite 
value corresponding to a localized configuration (super- 
conductivity). This is illustrated in Fig.|31for Ej/Ec = 1 
and a similar result is found for Ej/Eq = 2. In the sim- 
ulations, we used AtEq — 0.25, so the two coupling 
strengths correspond to AtEj = 0.25 and AtEj = 0.5 
respectively, and the inverse temperature is related to the 
lattice size by l3Ec = N/A. 

Another method to locate the phase boundary is to 
look at the Fourier transform {4i(j))ui„ of ((/)(O)0(r)) as was 
done in Refs. At low temperature, the Matsubara 

frequencies uin = are densely spaced and it becomes 
possible to determine the zero bias resistance R from an 
extrapolation to n = 0, 



R 1 

Hq n^O 2,71 



(9) 



The critical value of a can be determined from the tem- 
perature dependence of R. In the T = superconducting 
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FIG. 3: {{(f) — (j))'^) plotted as a function of system size (inverse 
temperature) for several values of dissipation strength a. The 
curves are parabolic fits to the data obtained for Ej/Ec ~ 1. 
The inset shows {{(f)— (f>)'^) as a function of a for Ej/Ec ~ 0.75 
and (from top to bottom) iV = 200, 80 and 36. 



phase, the resistance will decrease to zero with decreas- 
ing temperature. If the junction turns insulating, the 
resistance will increase and saturate at Rs, since normal 
electrons continue to flow through the shunt. In Fig. ^ 
we plot ■^\ujn\{(t>4>)i^„ for different temperatures and val- 
ues of a. At small n/N, the curves for a > 1.025 and 
their extrapolations to n = decrease with decreasing 
temperature, while those for a < 0.975 increase. 

In order to make the above analysis more precise and 
quantitative, we performed the extrapolation (jSJ using 
parabolic fits to the first five Matsubara points. The val- 
ues of R obtained by this procedure are plotted as a func- 
tion of the inverse temperature in Fig.El Again, the tran- 
sition at a = 1 (i? decreasing with inverse temperature 
above and increasing below the critical value) is clearly 
visible. Figs. 0] and |5l show the data for Ej/Ec ~ 1 but 
similar results were obtained for Ej/Ec — 2. This is the 
first convincing numerical evidence that the phase tran- 
sition occurs at a = 1.00(2) for intermediate coupling 
strengths, which supports the conjecture that ac does 
not depend on the value of the Josephson coupling. 

For a > 1, Korshunov |^ predicted a temperature de- 
pendence of the resistance of the form 



R{T) - 



(10) 



In a recent numerical study , 17 ] a temperature depen- 
dence incompatible with Eq. pO|l was found and the au- 
thors even observed a strong dependence of the power- 
law exponent on the Josephson energy. However, these 
data were obtained at rather high temperature, where 
Eq. (|10|l may not be valid. At the much lower tempera- 
tures accessible with the cluster algorithm, we find a good 
agreement with Eq. HlOfl . For Ej/Ec = 1 this is illus- 
trated in Fig. [SI by the curves proportional to iV~(^"~^^ 
A similar agreement was found for Ej/Ec = 2. 
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FIG. 4: ■^\uJn\{<j)4>)ui„ plotted as a function of n/N for dif- 
ferent system sizes (inverse temperatures) and dissipation 
strengths a. Ej/Ec = 1 and the values of a are (from top to 
bottom): 0.9, 0.95, 0.975, 1, 1.025, 1.05, 1.1 and 1.2. 
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FIG. 5: Zero bias resistance as a function of temperature for 
Ej/Ec = 1 and several values of the dissipation strength. 
From top to bottom, the data points correspond to a = 0.9, 
0.95, 0.975, 1, 1.025, 1.05, 1.1, 1.2, 1.3, 1.4, 1.5 and 1.6. The 
lines show a fit of the power- law Eq. lITOl . 



In this Letter we presented a new algorithm for the 
siniulation of resistively shunted Josephson junctions and 
similar physical systems. The algorithm combines local 
updates in frequency space with rejection-free cluster up- 
dates which exploit the symmetries of the cosine poten- 
tial. The cluster moves dramatically improve the sam- 
pling efficiency, allowing us to study junctions at tem- 
peratures which are 100 times lower than what has pre- 
viously been possible. This has enabled us to reliably 
confirm conjectured properties of the Schmid transition. 

This algorithm will be essential in the study of Joseph- 
son junction arrays. Already for the two-junction model 
with charge relaxation |l3| an interesting phase is pre- 
dicted to occur, which exhibits superconducting phase 
coherence between the leads, while the individual junc- 



tions are insulating. At intermediate values of Ej/Ec 
the quantum phase transition might be controlled by an 
intermediate coupling fixed point with as yet unknown 
properties. For extended systems, such as chains and 2D 
arrays, floating phases and lines of fixed points with con- 
tinuously varying exponents have been predicted 
Our algorithm will enable the numerical investigation of 
these models. 

We acknowledge support by the Swiss National Sci- 
ence Foundation and helpful discussions with H. G. Ev- 
ertz, Ch. Helm and G. Refael. This work was started at 
the Kavli Institute for Theoretical Physics, UCSB. The 
calculations have been performed on the Asgard Beowulf 
cluster at ETH Ziirich, using the ALPS library ji^. 
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